using dnAnalytics.Exceptions;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Solvers;
using NUnit.Framework;

namespace dnAnalytics.UnitTests.LinearAlgebra.Solvers
{
    [TestFixture]
    [Category("Managed")]
    public class BiConjugateGradientStabilizedSolverTest
    {
        private const double convergenceBoundary = 1e-10;
        private const int maximumIterations = 1000;

        [Test]
        [ExpectedException(typeof (NotSquareMatrixException))]
        public void SolveWideMatrix()
        {
            SparseMatrix matrix = (SparseMatrix) MatrixBuilder.CreateMatrix(2, 3, MatrixType.Sparse);
            BiConjugateGradientStabilizedSolver solver = new BiConjugateGradientStabilizedSolver(matrix);
        }

        [Test]
        [ExpectedException(typeof (NotSquareMatrixException))]
        public void SolveLongMatrix()
        {
            SparseMatrix matrix = (SparseMatrix) MatrixBuilder.CreateMatrix(3, 2, MatrixType.Sparse);
            BiConjugateGradientStabilizedSolver solver = new BiConjugateGradientStabilizedSolver(matrix);
        }
        
        [Test]
        public void SolveUnitMatrixAndBackMultiply()
        {
            // Create the matrix
            Matrix matrix = MatrixBuilder.CreateIdentityMatrix(100);
            
            // Create the y vector
            Vector y = VectorBuilder.CreateVector(matrix.Rows, 1, VectorType.Dense);

            // Create an iteration monitor which will keep track of iterative convergence
            ConvergenceMonitor monitor = new ConvergenceMonitor(maximumIterations,
                                                                convergenceBoundary,
                                                                matrix.InfinityNorm());
            BiConjugateGradientStabilizedSolver solver = new BiConjugateGradientStabilizedSolver(matrix, monitor);

            // Solve equation Ax = y
            Vector x = solver.Solve(y);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            Vector z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.AreEqual(SolutionStatus.Converged, monitor.SolutionStatus, "#04");

            // Now compare the vectors
            for (int i = 0; i < y.Count; i++)
            {
                Assert.AreEqual(y[i], z[i], convergenceBoundary, "#05-" + i.ToString());
            }
        }

        [Test]
        public void SolveScaledUnitMatrixAndBackMultiply()
        {
            // Create the matrix
            Matrix matrix = MatrixBuilder.CreateIdentityMatrix(100);
            // Scale it with a funny number
            matrix.Multiply(System.Math.PI);

            // Create the y vector
            Vector y = VectorBuilder.CreateVector(matrix.Rows, 1, VectorType.Dense);

            // Create an iteration monitor which will keep track of iterative convergence
            ConvergenceMonitor monitor = new ConvergenceMonitor(maximumIterations,
                                                                convergenceBoundary,
                                                                matrix.InfinityNorm());
            BiConjugateGradientStabilizedSolver solver = new BiConjugateGradientStabilizedSolver(matrix, monitor);

            // Solve equation Ax = y
            Vector x = solver.Solve(y);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            Vector z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.AreEqual(SolutionStatus.Converged, monitor.SolutionStatus, "#04");

            // Now compare the vectors
            for (int i = 0; i < y.Count; i++)
            {
                Assert.AreEqual(y[i], z[i], convergenceBoundary, "#05-" + i.ToString());
            }
        }
        
        [Test]
        public void SolvePoissonMatrixAndBackMultiply()
        {
            // Create the matrix
            SparseMatrix matrix = new SparseMatrix(100);
            // Assemble the matrix. We assume we're solving the Poisson equation
            // on a rectangular 10 x 10 grid
            int gridSize = 10;
            // The pattern is:
            // 0 .... 0 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 0 0 ... 0
            for (int i = 0; i < matrix.Rows; i++)
            {
                // Insert the first set of -1's
                if (i > (gridSize - 1))
                {
                    matrix[i, i - gridSize] = -1;
                }
                
                // Insert the second set of -1's
                if (i > 0)
                {
                    matrix[i, i - 1] = -1;
                }
                
                // Insert the centerline values
                matrix[i, i] = 4;
                
                // Insert the first trailing set of -1's
                if (i < matrix.Rows - 1)
                {
                    matrix[i, i + 1] = -1;
                }
                
                // Insert the second trailing set of -1's
                if (i < matrix.Rows - gridSize)
                {
                    matrix[i, i + gridSize] = -1;
                }
            }

            // Create the y vector
            Vector y = VectorBuilder.CreateVector(matrix.Rows, 1, VectorType.Dense);

            // Create an iteration monitor which will keep track of iterative convergence
            ConvergenceMonitor monitor = new ConvergenceMonitor(maximumIterations,
                                                                convergenceBoundary,
                                                                matrix.InfinityNorm());
            BiConjugateGradientStabilizedSolver solver = new BiConjugateGradientStabilizedSolver(matrix, monitor);

            // Solve equation Ax = y
            Vector x = solver.Solve(y);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            Vector z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.AreEqual(SolutionStatus.Converged, monitor.SolutionStatus, "#04");

            // Now compare the vectors
            for (int i = 0; i < y.Count; i++)
            {
                Assert.AreEqual(y[i], z[i], convergenceBoundary, "#05-" + i.ToString());
            }
        }
    }
}